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Abstract 

It is shown that the Mott insulating and superfluid phases of bosons in an optical lattice 

may be distinguished by a non-local 'parity order parameter' which is directly accessible 
via single site resolution imaging. In one dimension, the lattice Bose model is dual to 
a classical interface roughening problem. We use known exact results from the latter to 
prove that the parity order parameter exhibits long range order in the Mott insulating 
phase, consistent with recent experiments by Endres et al. [Science 334, 200 (2011)]. In 
two spatial dimensions, the parity order parameter can be expressed in terms of an equal 
time Wilson loop of a non-trivial U{1) gauge theory in 2 + 1 dimensions which exhibits 
a transition between a Coulomb and a confining phase. The negative logarithm of the 
parity order parameter obeys a perimeter law in the Mott insulator and is enhanced by 
a logarithmic factor in the superfluid. 
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1. Introduction 

A central assumption underlying both classical and quantum physics is that its basic 
laws are local, i.e., that the fundamental equations can be expressed in terms of relations 
between physical observables at a given point in space and time [1]. Symmetries and 
diS^erent types of order are thus related to the behavior of correlation functions of local 
observables. In particular, a qualitative change in macroscopic properties is typically 
associated with the appearance of long range order in some local observable 0{x). Its 
two-point correlation 

{dix)diy))^ 0^(00)^0 (1) 

thus approaches a finite constant as |a: — y| goes to infinity. In recent years, a lot of 
interest has focussed on systems where this standard characterization of different phases 
of matter by finite values of some local order parameter 0{x) fails. This is the case, e.g., 
in the Quantum Hall Effect, a paradigmatic example of a topological insulator. It is an 
incompressible state described by a Chern-Simons theory which has gapless excitations 
associated with chiral edge currents [2, 3]. 
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Our aim in the present paper is to study non-local orders that may be used to char- 
acterize Mott insulating phases of lattice bosons. Quite generally, Mott insulators are 
defined by their incompressibility n = dn/d(i = Q, & response function that is not as- 
sociated with long range order in any local observable 0{x). Specifically, we focus on 
Mott insulators which do not break any lattice symmetries due to, for instance, the 
formation of a commensurate charge density wave [4] or due to Neel order in the spin 
degrees of freedom, as happens in Mott-Heisenberg insulators [5] realized in undoped 
high-temperature superconductors [6] . To investigate possible non-local orders that may 
exist in otherwise featureless Mott insulators, the particular case of bosons in an optical 
lattice is of special interest. By a simple change of the lattice depth, they may be tuned 
to undergo a superfluid (SF) to Mott insulator (MI) transition [7]. Moreover, thanks 
to the direct accessibility of the atomic distribution by optical imaging methods which 
allow measuring in situ density distributions [8] and even arbitrary density correlations 
at the single- atom level [9, 10], they also provide a new perspective on the microscopic 
details of the involved states. The question of whether some non-local order exists in the 
MI phase of bosons in an optical lattice has been addressed before by Berg et al. [11] 
in the particular case of one dimension. While their main focus has been the study of 
hidden 'string-order' that appears in the presence of repulsive interactions of finite range, 
they have shown via Bosonization that even in the much simpler situation of pure on-site 
repulsion, there is a non-local observable which takes finite values in the MI and is zero 
in the SF. Using single site imaging, the associated 'parity order parameter' (POP) has 
been measured by Endres et al. [12]. For lengths of up to eight lattice spacings their 
results were consistent with the theoretical expectation of a non-vanishing parity order 
in the MI and an algebraic decay to zero in the superfluid phase. 

Our aim in the following is a detailed study of parity order in d = 1 and also in 
d = 2 dimensions using duality transformations. Specifically, in = 1, a number of 
exact results for the parity order parameter are obtained from a systematic expansion 
around the atomic limit and from analytical results for the roughening transition of the 
equivalent classical interface model in two dimensions [13, 14, 15]. In the d = 2 case, the 
lattice bosons are dual to a three-dimensional C/(l) lattice gauge theory. This mapping 
has first been discussed by Peskin [16] starting from the classical three-dimensional XY 
model and has been extended to the two-dimensional quantum XY model by Fisher 
and Lee [17]. As will be shown here, the parity order parameter is mapped in the dual 
model on a quantity which is 'more local by one dimension': in d = 1, it is mapped on 
a two-point correlation function for a local operator of the type in Eq. (1) which shows 
an algebraic decay in the SF phase and converges to a constant in the MI. In the d = 2 
case, the POP is mapped onto an equal time Wilson loop in the dual gauge theory, 
which serves to distinguish the superfluid and Mott insulating phases according to the 
dependence on the system size L. In particular, we flnd a perimeter law dependence in 
the MI phase which leads to an exponential decay of the POP while in the SF phase a 
logarithmic correction to the perimeter law leads to a super-exponential decay. 

2. Bose-Hubbard Model and Parity Order 

Ultracold bosons in an optical lattice provide a realization of the Bose-Hubbard model, 
as suggested theoretically by Jaksch et al. [18] and first realized experimentally by Greiner 
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et al. [7]. The associated many-body Hamiltonian [19] 

^BH = -JY1 + ^ XI ~ ^) ' 

{iJ) ' 

describes the coiiipetitioii between a kinetic energy which involves hopping to nearest 
neighbor lattice sites with amphtude J > and an on-site repulsive interaction U > 
which leads to an increase in energy if atoms hop to sites which are already occupied. 
Specifically, the operator a,; destroys a boson in a single particle Wannier state localized 
at lattice site i and the associated occupation number operator rij = d-aj has eigenvalues 

0, 1,2, We consider this model in both one and two spatial dimensions, specializing 

to the case of a simple quadratic lattice in the latter case. In either dimension, the ground 
state of the Bose Hubbard model exhibits a continuous transition between a superfluid 
and a Mott insulating state which may be tuned either by changing the ratio J/U at 
fixed density or by changing density at fixed J/U via the dimensionless chemical potential 
/u/f/ [19]. The universality class of this quantum phase transition is different in both 
cases. For the density driven Mott transition, the associated dynamical scaling exponent 
is z = 2 [20]. By contrast, changing the ratio J/U at fixed integer density n = 1,2,..., 
the transition at the tip of the Mott lobes has z = 1, i.e., it is described by an 0(2)-model 
in D = d + 1 dimensions which possesses a formal relativistic invariance [20]. It is this 
type of transition which will be considered in the following. 

The superfluid phase of the Bose Hubbard model is characterized by a conventional, 
local observable 0{x) = d^ associated with long range order in the one-particle density 
matrix, with O^(oo) = no as the condensate density^. This order can be observed in 
a direct manner in time-of-flight images [7], which measure the momentum distribution 
as the Fourier transform of the one particle density matrix [21]. The excellent quantita- 
tive agreement between the measured absorption images after time-of-flight and precise 
quantum Monte Carlo calculations which include both the effects of finite temperature 
and the harmonic trap potential [22], shows that the Bose- Hubbard Hamiltonian (2) 
provides a faithful description of cold atoms in an optical lattice. In the Mott insulating 
phase, the correlation function of the local bosonic field operator vanishes exponentially 
like (d^(x)d(O)) ^ exp(— |x|/^) with a correlation length ^ ^ 1/A that diverges like the 
inverse of the Mott gap A. The Mott phase is characterized by its incompressibility 
and has no associated order parameter, evolving in a continuous manner from a ther- 
mally disordered state as the temperature is lowered below the Mott gap A. Note that 
the observation of peaks in the noise correlations {fi{x)n{x')) in the Mott phase after 
time-of-flight by Foiling et al. [23] does not reflect any long range order: they are a con- 
sequence of the fact that the Fourier components YjR exp[i(fe — k') ■ R] of the average 
density tir = {fin) are of order N at wave vector differences k — k' = m{x — x')/ht 
which are equal to a reciprocal lattice vector G. 

As will be discussed in the following, the fundamental difference between the super- 
fluid and Mott insulating phase in terms of their compressibility implies that there is 



^This would not be true if the bosons were charged as, e.g., the Cooper pairs of a conventional 
superconductor, where the role of 0{x) is expected to be taken by the bi-Fermion operator tl^f^i {x). 
The associated correlation function {0'^(x)0{y)), however, is not gauge invariant and therefore does not 
constitute a proper order parameter, as noted by Wen [2]. 
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a non-local observable which behaves in a characteristically different manner in both 
phases. By a straightforward generalization of the parity order in Id introduced by Berg 
at al. [11], the non-local order is defined as 



where D is a spatial domain, i.e., an interval in d = 1 and an area in d = 2, and 

Pi = (—1)"* is the parity operator on lattice site i. There are two important proper- 
ties of this observable which should be noted right away: First of all, the observable is 
easily accessible in experiments since, due to light-induced collision losses, quantum gas 
microscopes directly measure the on-site parity rather than the actual occupation num- 
bers [10, 9]. As a second point, the observable must be calculated and measured in an 
open domain V which is part of a larger system, otherwise (O^('D)) = 1 would trivially 
be equal to one due to conservation of particle number. To study the dependence on the 
size, we shall characterize the domain V by its linear extension L measured in units of 
the lattice spacing. In analogy to the standard definition (1) of long range order, the 
parity order parameter (POP) {0'^{L)) exhibits long range order if (©^(oo)) is finite. 

This kind of order parameter is analogous to those studied in Ising models with a 
local gauge invariance, which have no conventional phase transitions to states with long 
range order, yet may exhibit different phases distinguished by non-local order parameters 
[24, 25, 26] (see also section 8). In the context of cold atoms, a more complicated 'string 
order' parameter was introduced, which characterizes a Haldane insulator that can form 
in one dimensional systems with longer range interactions [27]. Our focus is on the 
behavior of the parity order parameter at the conventional SF-MI transition, not only in 
Id [11, 12] but also in 2d. In particular, we will use a duality transformation to show that 
in two dimensions parity order is related to an equal time Wilson loop in a non-trivial 
U{1) gauge theory which exhibits different behavior as a function of system size L in the 
MI and SF phases. 

3. Number fluctuations and area law 

To obtain a qualitative understanding of the dependence of the non-local order defined 
in (3) on the size L of the domain we start by giving some qualitative arguments for the 
expected scaling behavior deep in the MI, where J/U <^ 1. Starting with the atomic 
limit, it is obvious that (0^(L))|j=o = 1 in any dimension since particle fluctuations 
axe then completely frozen. For small but finite J/U and for the simple case of a MI 
with average density n = 1, particle number fluctuations appear as pairs of empty and 
doubly occupied sites. Now, for an arbitrary number of pairs which are completely 
inside the domain V, the associated two minus signs in the product in Eq. (3) cancel. 
Only those pairs which are separated by the domain boundary lead to a reduction of 
(0^(L)) [12]. Intuitively one thus expects the parity order {0'^{L)) ~ exp(-L'^~^) to 
scale exponentially with the area of the domain's boundary. This intuitive expectation 
is supported by a systematic pcrturbativc calculation in an expansion around the atomic 
limit J/U = which — as will be shown in section 4 below — yields 
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Figure 1: (Color online) Illustration of how the order parameter is reduced from unity for d = 2. The 
grid lines indicate the lattice, bosons are represented by small circles. The domain 15, here taken to 
be a square, is shaded in gray. The minus signs from pairs which are completely inside or outside the 
domain (yellow ellipses) cancel out so that there is no contribution while pairs which are separated by 
the domain boundary (red ellipse) contribute a minus sign which leads to a reduction of (C'^(L)). 

up to second order in J/U <C 1. Clearly the expansion is well defined only in d = 1, 
while in higher dimensions the effective expansion parameter {J/U)'^ L"^^^ is small only 
up to system sizes of order (U/ jy^^'^~^K 

Further insight into the origin of the perimeter law for the decay of the parity order 
can be gained by assuming that the expectation value in (3) may be calculated within a 
Gaussian approximation such that 

where Sfii = — n. Within this approximation. — ln(0^(L)) is simply a measure of the 
total number fluctuations {SN'^) in a domain of size L as part of an infinite system. Now, 
the standard thermodynamic relation {SN"^) = ksT dN{ii)/d^i in this effectively grand 
canonical situation seems to indicate that these fluctuations vanish at zero temperature 
which would imply a trivial result (0^(L)) = 1 in the Gaussian approximation. This 
is not true, however, because the relation only applies in the thermodynamic limit and 
neglects boundary terms. For a careful calculation of {SN'^) at zero temperature and in a 
finite system, we generalize the analysis of Giorgini et al. ([28], see also [29, 30[) for Bose 
gases in a d = 3 continuum to arbitrary dimensions. The particle number fluctuations 

{SN^) =Sd drr'^-^T{r)ni^{r) , (6) 
Jo 

in a spherical domain of radius L can then be calculated from the pair distribution 
function 

i.(r) = <5(r)+n(.9(2)(r.)-l) = J ^e*'^-5(g) (7) 
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and the volume T(r) of the intersection between two d-dimensional balls of radius L 
separated by the distance r. Here, is the surface of a unit sphere in d dimensions 
while S{q) is the standard static structure factor. 

In the superfluid, the zero temperature static structure factor has the non-analytic 
behavior S{q) ~ Oi\q\ at small wave numbers, characteristic for any compressible phase. 
Since collective excitations exhaust the /sum rule for long wavelengths, the prefactor 
a = h/2mcs is, moreover, completely fixed by the exact sound velocity Cg. As a result of 
the non-analytic behavior of S{q), the associated pair distribution function exhibits an 
algebraic decay i^{r) ~ — a/vrr^ in d = 1 and fir) ~ —a/2TTr^ in = 2 at long distances. 
The behavior of the number fluctuations for large L is then found to be 

{5N^) ^ anL'^-Hn{L/^h) (SF) , (8) 

where the effective healing length = h/mcs serves as a short-distance cutoff. Within 
the Gaussian approximation, therefore, the parity order 

decays to zero with a power law in the superfluid phase in one dimension while in two 
dimensions, the decay is super-exponential. 

In the MI phase, the quadratic number fluctuations are much smaller and the scaling 
of {SN"^) with system size L differs in a qualitative manner from that in the SF. Indeed, at 
zero temperature, incompressibility of the MI implies that the structure factor vanishes 
analytically for 5 and thus S{q) = liq^"^ + ■ ■■ to leading order in an expansion 
around q = 0. Here, ^ is the characteristic length of the exponential decay of the 
one-particle density matrix mentioned in the previous section, while 7 is a numerical 
constant which relates the scales appearing in the density correlation and in the one- 
particle density matrix. From perturbation theory, one finds that in one dimension 7 is 
proportional to (J/U)^ to leading order [31]. 

Since the static structure factor is analytic around g = 0, the pair correlation function 
i/(r) decays exponentially on the characteristic scale ^, effectively cutting off the integral 
in Eq. (6). As a result, one finds that 

{5N^) hL"^-^ (MI) (10) 

scales like the area of the boundary of the domain with a coefficient h = bo'^ffi^^, where 
foo is a numerical constant depending on the geometry and the precise functional form of 

The results obtained within this Gaussian approximation are — of course — not exact. 
Yet, as will be shown in sections 6 and 7 below, they give the correct qualitative behavior 
of the parity order parameter as a function of the size L of the domain. On a basic 
level, therefore, the different behavior of the non-local parity order in the MI and the 
SF is a simple consequence of the fundamental difference between number fiuctuations 
in an incompressible versus a compressible system [32, 33[. In this context, it is also 
instructive to note that the underlying scaling (JTV^) ~ anL'^~^ ln{L/^) of the number 
fiuctuations in the compressible and gapless superfiuid compared to {5N^) ~ L'^"^ in the 
incompressible and gapped MI are reminiscent of similar results obtained for the scaling 
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of the entanglement entropy S{L). The fact that the fluctations of conserved quantities 
are closely related to the latter, exhibiting a simple area law for gapped phases, has 
been noted by Swingle and Senthil [34]. There are, however, a number of important 
differences: the entanglement entropy typically scales like S{L) = bL'^~^ with a non- 
universal prefactor b even in gapless phases while the corresponding number fluctuations 
(SN^) have an additional logarithmic enhancement factor for any phase with a finite 
compressibility. For the entanglement entropy, violations of the area law by a logarithmic 
factor of the form S{L) ~ L'^~^ lii(^) appear, e.g., in free fermions with a Fermi surface 
and also for Landau Fermi liquids, not in a gapless superfluid, however. Indeed, as noted 
by Metlitski and Grover [35], its entanglement entropy S(L) = bL'^^^ + AS* exhibits an 
additive — not multiplicative — logarithmic contribution /S.S = In (psL'^~^/cs)/2 which is 
universal, i.e., it only depends on the superfluid stiffness and the sound velocity Cg as 
effective low energy constants. 

The simplicity of the Gaussian approximation for studying the parity order also allows 
for a straightforward discussion of how the above results are affected by a finite temper- 
ature. For a neutral SF, the structure factor at temperatures ksT ^ mc^ reads [36] 
S{q) ~ Oi\q\ coth(cs|q|/2A;Br). The characteristic length scale at which the zero temper- 
ature result S{q) ~ a\q\ is equal to the q = Q thermal value S{q = 0) = ^k^Tajcs is 
therefore given by =■ hcs/2i:k^T = A|,/^h, where At is the thermal wavelength and 
— h/mCs the effective healing length. Since the parity order effectively probes number 
fluctuations at a finite wave vector q ~ 271 /L the zero temperature results remain valid 
as long as tt ^ L. Since typical values of the healing length are of order ~ 1 Hni, the 
characteristic scale of becomes larger than a lattice spacing at temperatures below 



An equivalent reasoning can be applied in the Mott insulating phase, where the 
thermal behavior of the structure factor S{q = 0) = nk-sTnT involves the compressiblity 
kt- For dimensional reasons, the latter must be of the form kt = (" A)^^/(/3A) with 
the Mott gap A and a function f{(3A) ~ e~^^ which has a thermally activated form. 
The zero temperature results are then valid as long as the domain size L satisfies 



In the interesting regime T <C A this is easily satisfied due to the exponential form of 
/(/3A). The actual constraint is rather the condition /3A ^ 1 which becomes increasingly 
hard to satisfy as one approaches the critical point. 

4. Perturbative analysis in the MI 

Deep in the MI, where J/U <C 1, exact results for the parity order can be obtained 
by a systematic perturbation theory around the atomic limit J = of the Bose-Hubbard 
model. Here we make use of a method developed by van Dongen [37], which is an 
extension of a formalism due to Harris and Lange [38]. The basic idea is to construct 
a canonical transformation of the creation and annihilation operators, = e^b^e~^ 
with an anti-hermitian generator S such that in the new basis, the occupation numbers 
and hence the interaction part remain invariant under hopping. Explicitly, one writes 



70 nK. 




(11) 
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the Hamiltonian in the form H = UD + K and defines the hopping term for the new 
particles 

f = - J ^ StSj- ^ k = e^fe-^ . (12) 

The requirement that D is invariant under liopping is obeyed provided it is a constant 
of motion, i.e., [-ff, -D] = 0. This condition fixes the transformation S, which may be 
calculated order by order in a systematic expansion 

^=Et^^"' (13) 

n>l 

in powers of 1/U . By expanding the exponential and substituting the expansion of the 
operator S, the POP takes the form 



{0'^{L)) = l + e-''^^"] 



($0 



k\ ^ °' 2! 



u 



C/2" 



u 



Si, 



$o) + • • • I , (14) 



where N{L) is the total number operator within the domain 2? and Nq is the eigenvalue 
of N{L) in the atomic limit ground state |<I>o)- Since |(f>o) is an eigenstate of N{L), the 
first term in the curly bracket vanishes at all orders. As a consequence, to calculate the 
POP at order n in J/U, one needs to determine the operators S'l, . . . , iSn-i- 

Specifically, substituting 5*1 from equation (A.l) (cf. Appendix A) and evaluating 
the commutators, one finds the result (4) stated already in section 3. In one dimension, 
we have continued the perturbative expansion up to fourth order in J/U. For reasons of 
space, the details of this calculation are deferred to Appendix A and we only give the 
result here: 



(02(L)) = l-8n(n + l) 



- -n(n+ l)[n(473n 



217) - 234] 



+ ■ 



(15) 



Just as the leading order term, the next-to-leading correction is independent of the 
domain size L as long as the latter is larger than the order of pertTirbation. Hence, 
{0'^{L)) is independent of L for sufficiently small J/U, in agreement with a DMRG 
calculation in [12] where {0'^{L)) was found to be essentially independent of L for J/U < 
0.1 for domain sizes ranging from 1 to 60. 



5. Duality transformation 

In the following, we show that within a reduced description of the SF to MI transition 
at fixed density in terms of a quantum rotor model [20], exact results for the parity order 
parameter can be obtained from analyzing a (d + l)-dimensional classical lattice model. 
This is particularly interesting in the d = 2 case where the parity order {0^{L)) can be 
expressed in terms of an equal time Wilson-loop in a nontrivial U{1) gauge field which 
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is dual to the original lattice boson model. Our aim is twofold: first of all, unlike the 
qualitative arguments presented in section 3, the duality transformation permits us to 
derive exact results for the large L behavior of in both the MI and SF phases 

within one unified framework. From a different point of view, however, the duality of the 
lattice boson model to a dynamical gauge field in 2 + 1 dimensions can also be read the 
other way, where ultracold atoms in an optical lattice provide a quantum simulator of a 
nontrivial lattice gauge theory. 

The starting point for the mapping is based on the realization that the SF to MI 
transition at fixed density is driven by phase fluctuations only [19]. To isolate these from 
the full Bose-Hubbard Hamiltonian, it is convenient to consider the limit of large filling 
n 3> 1, where the boson operators may be rewritten in a density-phase representation, 
cij ~ e^'^i i/n + Sfij [19]. In the limit n 1, the number fluctuations 5hj can be expanded 
up to second order and the Bose-Hubbard Hamiltonian is thus transformed into the 
Hamiltonian 

Hj = ^Y.^l + EiY,[^- cos((^,+„ ~ ^,)] (16) 

X x,u 

of a system of quantum rotors at each lattice site x with a discrete eigenvalue spectrum 
0, ±1, ±2, . . . which are coupled to nearest neighbors a; + u by a Josephson energy Ej = 
2nJ. (Note that we have redeflned 5hx — n^, i.e., is now the deviation from the 
average local occupation number). 

Following a standard procedure, the partition function associated with the quantum 
rotor model (16) can now be expressed in terms of a path integral by dividing the interval 
[0,/3] into A'^^ subintervals of width e = P/Nr- Inserting a complete set of number states 
at each step, one thus obtains a discretized path integral representation of the form 

Z,= J2 i{n=^A}\e-'"'\{nx,2})---{{n^,NA\e-^"^\{n.A})^ (17) 

where the notion {n^.j} is a reminder of the fact that — at given j = 1, . . . ,iV^ — there 
are N'^ variables Uxj e Z for a; e (1, . . . , AT)** that have to be summed over all integers 
Z. In the limit £ — )• 0, the non-commuting terms in exp{—eHj) can be factorized to give 

The matrix elements of the Josephson coupling terms can now be simplifled by using the 
so called Villain approximation^ [39] 

exp{-££'j[l - cos{(f)x+u - 0x)]} - exp 1 - imx,u{4>x+u - <i>x)] ■ (19) 

Since the fundamental symmetry — > (j) + 2n due to the discreteness of the boson number 
is retained, this leaves the physics qualitatively unchanged. In one spatial dimension, this 



constant prefactor in (19) is suppressed because it only gives an irrelevant overall shift of the free 
energy. 
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approximation introduces one integer niij per lattice site, in two spatial dimensions, it 
introduces two integers trixj = {m^.^^j^mx^yj) at each lattice site. 
For each given j, one now uses the fact that 



JJ exp {—im, 



'x,u,j \^^x-\-u 



m 



X — U,U,J } Yx 



(20) 



for a periodic chain with 4'x+Nu — 4'x and the identity (7i'|e'™'^|Ti) — ^n' ,n+m since the 
operator e""''* shifts the particle number by m. The resulting partition function 



Zj = ^ exp 



eU \ - 2 

-^l^^x^j 

x,j 



2eEj ^ ""^'"'J' 

X,U,j 



n ^ 



(21) 



has a Gaussian form, however the variables nx,j and rrixj are integer-valued and are 
connected by the constraint Vaj-m-t-Vi-n = 0, where V^^r denotes the discrete derivative 
on the dual lattice of links along the physical directions x and the 'time' direction r. 
Thus, the variables n^,] and rrix.j together form a divergenceless (d + l)-dimensional 
integer vector field n = (n, m). In d = 1, this constraint may be resolved by introducing 
a single integer field hx,T such that n = Vxh and m = —Vrh. The partition function 
then becomes that of the discrete Gaussian model with height variable h which describes 
the roughening transition of a 2d interface. Moreover, the POP is mapped on a two-point 
correlation function of the local variable 0{x) = exp(i7r/i(a;)). This model is discussed in 
section 6. 

In d = 2, the constraint is automatically satisfied if one introduces a three-component 
vector potential a such that n = V A a, where the lattice curl is defined as (V A ax)i = 

k^ijkio-x-k k ^ ^x-j-k fe)- "r^*^ partition function then becomes that of a (2-1- 1)- 
dimensional U{1) gauge theory, as will be discussed in detail in section 7. Remarkably, 
under the duality transformation the parity order parameter is mapped onto an equal 
time Wilson loop [40] 



{0\L)) = exp 



i-K ^ d^x (V A a)^ 
xev 



= { exp 



xedv 



(22) 



where we have used the discrete version of Stokes's theorem on the last line, dT> is the 
boundary of the domain V and Ax is a unit vector directed along the boundary in the 
positive mathematical sense. The Wilson loop is a gauge-invariant quantity which is 
often used to characterize phases in gauge theories [41]. As will be shown in section 7, 
in the present case, the transition between a phase with massless and one with massive 
'photons' in the underlying ?7(l)-gauge theory predicts qualitatively different behavior of 
the parity order in the SF and MI phases of the original lattice Boson model, consistent 
with the qualitative considerations discussed in section 3. 



6. Discrete Gaussian interface model 

As shown in the previous section, in the d = 1 case the partition function can be repre- 
sented in terms of a single integer h on each lattice site such that (n, m) = {Vxh, —Vrh). 
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Choosing e = the resulting partition function 



exp 



{hi,} 



1,3 ^ ^ 



(23) 



defines an anisotropic discrete Gaussian (DG) model for an integer valued height variable 
hij above a two-dimensional, perfectly flat interface hij = (note that an overall shift 
hij — > hjj + Z of this reference plane is irrelevant) . This mapping has been used earlier in 
the context of the SF-MI transition in one dimension by one of the present authors [42]. 

The DG model is a classical model which exhibits a phase transition from a smooth 
interface in the regime where U dominates to a rough phase in the limit E;] ^ U . The 
smooth phase, which corresponds to the MI in the original quantum rotor model, is 
characterized by a finite dimensionless step free energy fg [43] which is related to the 
Mott gap A/i of the dual model (16) by 2/s = A/i/C/ [42]. The dimensionless step free 
energy is a decreasing function of Ej/U and reaches /g = 1/2 at Ej = 0. In this limit 
an additional boson is described by a step of unit height which is parallel to the r axis, 
i.e., the boson world lines exhibit no quantum fluctuations. 

When Ej/U ~ 1, the model is essentially isotropic. To render this manifest, it is 
convenient to choose s = 1 /\/EjU so that 



exp 

{hij} 



1,3 



(24) 



In this representation, the ratio Tdg = \/Ej/U of kinetic and interaction energy in 
the underlying quantum rotor Hamiltonian (16) plays the role of an effective tempera- 
ture. The discrete Gaussian model (24) is known to have a roughening transition of the 
Kosterlitz-Thouless type [44] at a critical temperature Tr ~ 0.73. In the smooth phase, 
the mean square surface displacement 



^h\L) = {{hi^-hi,rf) 



(25) 



remains finite as the distance L = — between two points on the surface 

approaches infinity. By contrast, the rough phase of the discrete Gaussian model at 
Tdg > is characterized by a logarithmically divergent Ah^{L) ~ ln(L). 

These results on the discrete Gaussian model, for which the existence of a Kosterlitz- 
Thouless transition has been proven rigorously by Frohlich and Spencer [13], can now 
be translated back to understand the nature of non-local order in the original Bose- 
Hubbard or quantum rotor model. In particular, the qualitative results that were derived 
in section 3 within a Gaussian approximation from considering the number fiuctuations 
within a domain of linear size L can now be put on a rigorous footing. This relies on the 
fact that number fluctuations in the original model of bosons hopping on a lattice are 
transformed, via the duality, to fluctuations of the normal vector of the 2d interface by 



^Note that the Trotter decomposition in Eq. (17) is usually performed for finite /3 and thus Nr = 
GO requires e 0. Here we keep e finite but consider the Umit of zero temperature /3 oo. 
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n = Vxh. As a result, the nonlocal order parameter defined in Eq. (3) translates into 
the characteristic function 

{0^{L)) = {eMi^m - h{Q)])) (26) 

of the probability distribution p{h, \ {l,j) — {I' = {^{hij — hfj' — h)) that the height 
variables at two sites at a distance L differ by h. 

For a detailed understanding of the behavior of the parity order parameter, we can 
now use exact results on the classical roughening transition of two-dimensional interfaces. 
Specifically, within the smooth phase, Forrester [14] has obtained the exact probability 
distribution for the height of a lattice site near the center in the body-centered solid-on 
solid (BCSOS) model for fixed boundary conditions and the thermodynamic limit: it 
turns out to be a discrete Gaussian distribution p{h) = Afe~^'^~^^^^ 1'^'^ , where AA is a 
normalization constant, and with = C j — T/Tji where C = ^J2/h\2. The non- 
zero expectation of this distribution stems from the fact that in the BCSOS model, one 
has to deal with two sublattices where the sites of one take only even values while the 
sites of the other take only odd values. Since the outermost sites are fixed at values and 
1 (according to the sublattice), the height in the center may equivalently be seen as the 
height difference to the border, or, in the thermodynamic limit, as the infinite-distance 
limit of this difference. Hence, the characteristic function of this probability distribution 
is just the two-point correlation function discussed above within the DG model. Since 
the DG and the BCSOS model are in the same universality class, we may conclude that 
in the limit of infinite distance, the probability distribution for the height difference 
between two points in the DG model equally becomes a discrete Gaussian distribution 
p(Aft = n) = 7Ve-"'/2<T2 ^.^^^^ ^ Cdg/\/i - ^bc/Tk- Using the Poisson formula, 
one then obtains the characteristic function 

j:Pin)e"'- = ^ E e-^-^'^-^^-r/^ , (27) 

i.e., the characteristic function is a periodic sum of Gaussians centered on 0, ±27r, ±47r, . . . . 
The prefactor is just the ratio of the normalization constants of the continuous and the 
discrete Gaussian distribution and tends to unity for cr^ 1. As Tdg approaches Tr 
from below, diverges and the individual peaks of the characteristic function become 
very narrow. There are then only two contributions to the POP (from the peaks centered 
on and 27r) which goes to zero as 

(02(oo))=p(/c = 7r)~2exp I ^!^^ | , (28) 

V 2^1-fuG/fRj 

confirming the conjecture (0^(L)) ^ exp{— A[( J/C/)c — {J/U)]~^^'^} made in [12] on how 
the POP reaches zero as one approaches the critical point from the smooth phase (close 



*Note that a priori, Cdg 7^ C [14]. The numerical values of the roughening temperatures are model- 
specific as well. 
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to the critical point, writing {0'^{L)) as a function of J/U rather than Tdg oc \fjJU 
only affects the non-universal constant A). 

To understand the behavior of the parity order within the superfluid, which corre- 
sponds to the rough phase of the associated 2d classical interface model, it is convenient 
to use the equivalence between the DG model and the classical two-dimensional Coulomb 
gas (CG) derived by Chui and Weeks [44] . Within the 2d CG picture, the underlying SF- 
MI transition is translated into a phase transition of the Kosterlitz-Thoulcss type between 
an insulating phase where charges are bound — the SF phase of the original model — and 
a metaUic phase of effectively free charges which describes the MI. The metallic phase 
is characterized by a divergent polarizability J (Pr r'^p{r), where p(r) is the probability 
distribution function for the distance between a pair of opposite unit charges added to 
the system. Indeed, as shown by Chui and Weeks [44], for < ^ < 27r the two-point 
correlation function of the DG model maps on the partition function of the Coulomb gas 
in the presence of two opposite charges ±^ at a distance r = \ri — r2\, 

where F{r, ^) is the free energy of the neutral Coulomb gas in the presence of the added 
charges. It is well known [45] that when these added charges have unit strength [within 
the mapping by Chui and Weeks, a unit charge corresponds to ^ = 27r in F{r,^)], 

g-/3cGF(r,2,r) = ^ In = r"^'/^" , (30) 

with the dielectric constant £q which characterizes the insulating phase of the Coulomb 
gas. At the unbinding transition, the polarizability diverges, i.e., the critical value of the 

dielectric constant obeys (e^/co)c = 4. In our case of ^ = tt, we are dealing with a pair 
of half-unit charges^. As a result the POP 

{O^iL)) ~ r-''/^'° , (31) 

decays algebraically in the SF phase of the original lattice Bose model, and the mapping 
to the Coulomb gas gives us the exact value of the universal jump of the exponent at the 
roughening temperature: for Tdg — ^ from above, the exponent e^/4eo approaches one 
before jumping to zero, a consequence of the universal jump of the superfluid stiffness 
at the SF-MI transition of the Id quantum rotor model (for a more detailed discussion 
see [46]). 

7. U{1) gauge theory in 2 + 1 dimensions 

In d = 2, the duality mapping discussed in section 5 and the introduction of the three- 
component vector potential a with n = V A a lead to a (2 + l)-dimensional U(l) gauge 
theory [2, 16, 17, 47, 48] (for an alternative approach to this type of duality mapping 



^To avoid possible confusion, we stress that the right-hand side of Eq. (29) is well-defined for arbitrary 
^, but the equality to the left-hand side (i.e., the existence of the mapping) is restricted to values < 
^ < 27r, and thus to fractional charges. For J = 27rn with integer n, one trivially has e'^[''(''i^~'''''2^1 = 1. 
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(29) 



using the operator formalism, cf. [49]). Similarly to the Id case, it is convenient to choose 
£ = 1/ y/UE], which results in a partition function of the form 



^ = E-pM/5E(VAa)^) . (32) 

{a} \ ^ x,r / 

Since (V A a)^ = + 6^ is the energy density associated with a two-component 'electric 
field' Ci = Vrfli — Vittr and a scalar 'magnetic field' b = V^dy — Vya^, this partition 
function appears to describe the free field theory of pure electrodynamics in (2 + 1) di- 
mensions [2]. It is gauge invariant since the action only depends on the curl of a and is 
thus unchanged if the lattice gradient of an arbitrary function of x and r is added to a. 
What makes the model in Eq. (32) nontrivial is the fact that all fields take only integer 
values on a discrete space-time lattice. The SF to MI transition of the underlying lattice 
Bose model shows up at the level of the dual J7(l) gauge theory as a transition between 
a phase in which a really is a free field and the photon is massless and one with massive 
photons. 

In order to understand the physical meaning of this transition in terms of the gauge 
field degrees of freedom, it is convenient to consider the equal-time one-body density 
matrix in the quantum rotor model (16), which is mapped onto a ratio of two gauge field 
partition functions 



Mm = ^,ih^),-im^ = M = exp(-Ai=^[x,0]) , 



(33) 



where Z\x^ 0] differs from Z in that the constraint V • 'n-(y) = is replaced by V • n(y) = 
5y,x — Syfi [50]. Physically, this corresponds to a configuration with a pair of oppositely 
charged magnetic monopoles situated at x and 0. The fact that the one-body density 
matrix approaches a constant in the SF and decays exponentially in the MI thus leads 
to a fundamentally different behavior of the dimensionless free energy increase 

r , (MI) , , 

^ ^ [const. - 3^ (SF) ^ ' 

associated with the introduction of a monopole-antimonopole pair at distance x in the 
dual gauge theory (here, and Cg denote the superfiuid stiffness and sound velocity, 
respectively). In particular, the exponential decay of the one-body density matrix in the 
MI phase translates to a linear confinement while in the SF phase the monopoles interact 
via a 3d Coulomb potential. 

An effective low energy description of the gauge theory which properly accounts for 
the two different phases is provided by a Gaussian model of the form 



S=XJd'x!^[VAa{x)]' + ^[a{xrj , (35) 

where a is now treated as a continuous variable and T ~ Ej /U is a renormalized 
dimensionless temperature or coupling constant, similar to the one in the previous section. 
In the SF regime, where T is above a critical value of order one, the gauge field is in 
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its Coulomb phase where ^ = 00. Elementary excitations are then massless photons, 
which are just the phonons of the superfluid with linear dispersion to = Cs\q\ (note 
that 'photons' in (2 + l)-dimensional electrodynamics have no polarization degrees of 
freedom). By contrast, in the MI for small values of T, the gauge field is in a confining 
phase, with ^ = 1/m finite. The photons thus acquire a mass m and now represent the 
elementary particle-hole excitations of the MI with dispersion to = Cssjrn^c^ + (p (at 
the transition, one has Cg = 4.8J for the Bose-Hubbard model with n =- \ [51]). When 
evaluating expectations for the massless phase, it is convenient to keep a finite ^ during 
the calculation and only take the limit ^ — >■ 00 at the end of the calculation, which avoids 
the necessity of introducing an explicit gauge-fixing term [52] . This can be seen explicitly 
in the correlation function of the vector potential, which reads 



(a(q),a(q')fe)=r(27r)3(5(g + q') 



life 



+ 



(36) 



where \Ptisif\jk = ^jk — (lj(lk/<i^ and [Pi{q)]jk = QjQk/o^ are the components of the 
transverse and longitudinal projector with respect to q, respectively. Since the effective 
model (35) is Gaussian, the calculation of the parity order parameter from the Wilson 
loop in Eq. (22) is now easy and gives^ 



{0\L)) = exp 




/ d^a;[VAa 
J-D 




Using (36), the expectation appearing in the exponent in Eq. (37) reduces to 




d'^x[V Aa{x)]r 



T 



(2^ q2 + ^- 



(37) 



(38) 



where qj^ = {qx,Qy)- Note that q has three finite components whereas x is restricted to 
the T = plane. For a circular disk of radius R, the spatial integral appearing in (38) 
gives 



X 



2tvR 



Ji(5±i?) . 



(39) 



In the deconfincd phase, which corresponds to the superfluid, the mass parameter 1/^ is 
equal to zero and the resulting parity order parameter is given by^ 



\n{0^{R)) - 7r3fi?ln(7ri?) 



(40) 



in agreement with the qualitative result obtained in Eq. (9) within the Gaiissian approx- 
imation. Conversely, in the confined phase corresponding to the MI, The q^ is negligible 
compared with 1/^^, leading to a perimeter law 



-\n{0^{R)) r^f^^2TrR\ 



R 



(41) 



®The resulting scaling of the parity order will thus turn out to be the same as found in section 3 up 
to numerical factors, since both are obtained from a Gaussian model. Note, however, that the mapping 
of to a Wilson loop in the dual gauge theory provides a proper justification for these scalings. 

^The momentum integrals in (38) are cut off at A ~ tt since the original problem is defined on a 
lattice with lattice constant set equal to unity. 
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Thus, the characteristic decay length of the POP scales as ^ A^^U/J close to the 
critical point. Remarkably, a measurement of the parity order, which only involves the 
statistics of number fluctuations, therefore allows to extract the Mott gap. 

8. Conclusion &c Outlook 

We have presented a detailed study of parity order {0'^{L)) for lattice bosons in 
both one and two spatial dimensions, using duality transformations. Consistent with 
previous theoretical work [11] and recent experiments [12], it has been shown that the 
Mott insulating phase in one dimension exhibits long range parity order. An intuitive 
understanding of this result relies on the observation that for any incompressible phase, 
the number fluctuations {SN"^) ~ L'^~^ in a domain of size L scale with the area of the 
boundary ^ L'^^^. Using the duality to a discrete, classical interface roughening problem 
in two dimensions, these results have been put on a rigorous footing. In two dimensions, 
the parity order again shows qualitatively distinct behavior in the MI and SF phase. 
In this case, the variable (©^(i)) can be expressed in terms of an equal time Wilson 
loop of a nontrivial U{1) gauge theory in 2 + 1 dimensions. This is related to the fact 
that the density fluctuations in the original lattice Bose model are mapped to the scalar 
magnetic field in the dual gauge theory. A quite interesting result obtained from this 
mapping is the fact that the decay of parity order in the MI allows to measure the Mott 
gap from the statistics of number fluctuations. Since experimental measurements of the 
parity order in two dimensions are straightforward in principle, this seems a promising 
route to infer microscopic information from single site resolution imaging which is very 
diflicult to obtain otherwise. The detailed numerical factors in the scaling of {0^{L)) 
can unfortunately not be predicted from our effective long wavelength description of the 
quantum rotor model. Since we are dealing with a bosonic system, however, effective 
numerical methods are available to obtain precise results in a realistic experimental 
setup, as has been shown for thermodynamic properties and excitation energies [51]. 
In particular, numerical simulations directly deal with the Bose Hubbard model which 
applies in the relevant case of low flUing n = 1 instead of the qualitatively similar situation 
n ^ 1 that has been studied here within the quantum rotor model. 

An important open question is, of course, to which extent our results for Bose Mott 
insulators can be generalized to the fermionic case, which has also been realized experi- 
mentally with ultracold atoms [53, 54]. Based on the qualitative description in terms of 
the scaling of number fluctuations in section 3, we expect that the results obtained here 
carry over to the fermionic case despite the fact that no duality transformations exist in 
this case which allow to connect the parity order with a Wilson loop in a dual gauge the- 
ory. The presence of long range parity order for Id Mott insulators also in the fermionic 
case is consistent with the results of a recent DMRG study of the fermionic Hubbard 
model by Montorsi and Roncaglia [55]. A different kind of non-local order characterized 
by sub-lattice parity was found in this Model by Kruis et al. several years earlier [56, 57]. 

Finally, a quite interesting direction of further research on non-local orders for cold 
atoms in optical lattices is connected with the recent realization of a Id transverse Ising 
model using a tilted optical lattice [58]. Extending this setup to the case of two dimen- 
sions, a number of complex phases may appear depending on the type of lattice and the 
direction of the tilt [59]. A quite intriguing perspective would appear in a setup that 
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allows to realize a ferromagnetic version of this standard model for quantum phase transi- 
tions [20]. In the ferromagnetic case, the transverse Ising model is self-dual in one dimen- 
sion, while in two dimensions the dual theory is given by an Ising gauge theory [41]. For 
the latter, one can define non-local correlation functions as C(I?) = (Ilme-D ^m)^ where 
is the X component of the spin operator at site m and V is an area in two dimensions, 
in complete analogy with our definition of the parity order in equation (3). After the 
duality transformation, this observable is transformed into C(2?) = (Wk(^i-w ^k) • where 
(7^ is the z component of the spin operator at site k in the dual theory and dV are 
the two sites at both ends of the string in one dimension or the border of the area in 
two dimensions. As a result, one has in one dimension that limL^coC(i) > for the 
paramagnetic phase. In two dimensions, C(V) is a Wilson- Wegner loop around a closed 
path [24, 41], which shows an exponential scaling with the perimeter of the loop in the 
paramagnetic phase and with the enclosed area in the ferromagnetic phase of the orig- 
inal model [41, 60]. However, the detection of the non-local order parameter requires 
the measurement of the x component of the spin operator, which in turn requires the 
single-site resolved detection of the phase coherence between superposition states with 
different on-site occupation numbers, a technique that is so far not available. 
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Appendix A. Details about the perturbative calculation 

According to the program outlined in Section 4, calculting the POP to nth order 
in J/U amounts to calculating the operators S\,. . . ,Sn-i- Hence, for the second order 
result, we only need to calculate one operator. One finds 

= E '^^Po^tPo, , (A.1) 

DuD2 ^ ^ 

where the sums go over the eigenvalues D of the operator D and Pd is the projector on 
the eigenspace corresponding to this eigenvalue. Moreover, here and in the following, a 
prime on a sum indicates that values which make the denominator vanish are excluded 
from the sum. 

One may imagine the calculation of the expectation {0^{L)) as summing over all pos- 
sible hopping processes where the order in J/U indicates the number of occurring hops. 
For example, at second order, starting from a situation with uniform filling corresponding 
to |$o)) all possible processes consist of a single particle hopping to a neighboring site 
and back again. Different contributions stem from the position of the starting lattice site 
and its neighbor relative to the domain boundary and from the position of the operator 
iV'(D) in the product of operators. 
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The next non- vanishing contribution to the POP is of fourth order, so we need to 
calculate 82,3- The former turns out to be given by 



-Di,-D2,-D3 ^ \ ^ J J 1/ 

+ E \,2 (PD,fPn,fPn,-pD,fPn,fPn,) . (A.2) 

The expression for S3 is too unwieldy to reproduce here, so we Umit ourselves to a 
description of the involved hopping processes. They can be grouped into four types: 

• Back- and forth hoppings of two particles separated by more than two lattice sites, 
i.e., independent second order processes. 

• Processes where the same particle or hole hops twice before returning to its original 
site. 

• Processes where a particle hops twice and is then followed by the created hole or 
vice versa. 

• Processes where two particles start or end on the same lattice site before hopping 
back (cf. Fig. A.2). 



••^•^•9 9999^99 

999 999 ••••••• •••• 

••••••• ••••••• ••••••• 

••••••• ••••••• ••••••• 

Figure A.2: (Color online) Possible processes at fourth order where two particles start or end on the 
same lattice site. The images show the configuration after two hopping events (the remaining two restore 
uniform filling) for n = 4. 

Summing over all contributions finally yields Eq. (15). 
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